Please focus on association between factor 1 and striatum. I include calculations on data with and without imputation. Im glad to see they provide the same output, namely association between (subtype of) MSA and factor 1/striatum.

Setup

system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(magrittr)
library(dplyr)
library(MOFA2)
library(ggplot2)
library(UpSetR)
library(scHelper)
library(qs)
library(tidyr)
library(purrr)
library(tibble)
library(MOFAdata)

tt <- Sys.time()

Helper functions

Weights with information about imputation

plot_top_weights_imp <- function (object, imp.list, view = 1, factors = 1, nfeatures = 10, abs = TRUE, scale = TRUE, sign = "all") 
{
  if (!is(object, "MOFA")) 
    stop("'object' has to be an instance of MOFA")
  if (nfeatures <= 0) 
    stop("'nfeatures' has to be greater than 0")
  if (sign == "all") {
    abs <- TRUE
  }
  if (is.numeric(view)) 
    view <- views_names(object)[view]
  stopifnot(view %in% views_names(object))
  view <- MOFA2:::.check_and_get_views(object, view)
  factors <- MOFA2:::.check_and_get_factors(object, factors)
  W <- get_weights(object, factors = factors, views = view, 
                   as.data.frame = TRUE)
  if (scale) 
    W$value <- W$value/max(abs(W$value))
  W <- W[W$value != 0, ]
  W$sign <- ifelse(W$value > 0, "+", "-")
  if (sign == "positive") {
    W <- W[W$value > 0, ]
  }
  else if (sign == "negative") {
    W <- W[W$value < 0, ]
  }
  if (abs) 
    W$value <- abs(W$value)
  W <- W[with(W, order(-abs(value))), ]
  W <- as.data.frame(top_n(group_by(W, factor), n = nfeatures, 
                           wt = value))
  W$feature_id <- W$feature
  if ((length(unique(W$view)) > 1) && (nfeatures > 0) && (any(duplicated(W[W$factor == factors[1], ]$feature_id)))) {
    message("Duplicated feature names across views, we will add the view name as a prefix")
    W$feature_id <- paste(W$view, W$feature, sep = "_")
  }
  W$feature_id <- factor(W$feature_id, levels = rev(unique(W$feature_id)))
  
  # Add information on imputation
  W %<>% 
    mutate(feature_stripped = as.character(feature_id) %>% 
             strsplit("_") %>% 
             sget(1))
  
  W$imp <- imp.list[[view]] %>%
    .[match(W$feature_stripped, names(.))] %>% 
    factor(levels = c(F, T), labels = c("Not imputed", "Imputed"))
  
  se.assay <- SummarizedExperiment::assay(se.noimp) %>% 
    .[W$feature_stripped, grepl(view, colnames(.))]
  
  assay.per.cond <- se.assay %>% 
    as.data.frame()
  
  nsamples <- ncol(se.assay)
  
  # Add per condition
  assay_long <- se.assay %>% 
    as.data.frame() %>% 
    as_tibble(rownames = "protein") %>%
    pivot_longer(-protein, names_to = "id", values_to = "value") %>%
    left_join(meta %>% dplyr::select(id, Diagnosis), by = "id")
  
  diag_counts <- assay_long %>%
    group_by(protein, Diagnosis) %>%
    summarise(
      impn   = sum(is.na(value)),
      nsamp  = n(),
      impper = round(impn / nsamp * 100, 1),
      .groups = "drop"
    )
  
  diag_wide <- diag_counts %>%
    dplyr::select(protein, Diagnosis, impn, impper) %>%
    pivot_wider(
      id_cols    = protein,
      names_from = Diagnosis,
      values_from = c(impn, impper),
      names_sep  = "_"
    ) %>% 
    .[match(rownames(se.assay), .$protein),] %>% 
    dplyr::select(-protein)
  
  nsamples <- ncol(se.assay)
  
  W %<>%
    mutate(impn = se.assay %>% 
             apply(1, \(x) sum(is.na(x))) %>% 
             unname()) %>% 
    mutate(impper = round(impn / nsamples * 100, 1) %>% 
             unname()) %>% 
    cbind(diag_wide) %>% 
    dplyr::select(-feature, -view, -feature_stripped) %>% 
    mutate(feature_id = as.character(feature_id),
           imp = as.character(imp),
           factor = as.character(factor),
           impper = as.character(impper)) %>% 
    {rbind(data.frame(factor = "Factor1", value = 0, sign = "-", feature_id = "Columns", imp = "Not imputed", impn = "Tot, n", impper = "Tot%", impn_CTRL = "CTRL, n", impn_MSA = "MSA, n", impn_PD = "PD, n", impn_PSP = "PSP, n", impper_CTRL = "CTRL%", impper_MSA = "MSA%", impper_PD = "PD%", impper_PSP = "PSP%"), .)} %>% 
    mutate(feature_id = factor(feature_id, levels = rev(unique(feature_id))))
  
  if (any(is.na(imp.logi))) warning("There are NAs in the list of imputed proteins.")
  
  # Plot
  p <- ggplot(W, aes(x = feature_id, y = value, col = imp)) + 
    geom_point(size = 2) + 
    geom_segment(aes(xend = feature_id), linewidth = 0.75, yend = 0, show.legend = F) +
    scale_discrete_manual(aesthetics = "col", values = c("Imputed" = "firebrick", "Not imputed" = "black")) +
    coord_flip() + 
    labs(y = "Weight") + 
    theme_bw() + 
    theme(axis.title.x = element_text(color = "black"), 
          axis.title.y = element_blank(), 
          axis.text.y = element_text(size = rel(1.1), hjust = 1, color = "black"), 
          axis.text.x = element_text(color = "black"), 
          axis.ticks.y = element_blank(), 
          axis.ticks.x = element_line(), 
          legend.position = "right", 
          legend.title = element_blank(), 
          legend.text = element_text(color = "black"), 
          legend.key = element_rect(fill = "transparent"), 
          strip.text = element_text(size = rel(1.2)), 
          panel.background = element_blank(), 
          panel.spacing = unit(1, "lines"), 
          panel.grid.major.y = element_blank(), 
    ) #+ 
    # facet_wrap(~factor, nrow = 1, scales = "free")
  if (sign == "negative") 
    p <- p + scale_x_discrete(position = "top")
  if (abs) {
    p <- p + 
      ylim(0, max(W$value) + 0.1) + 
      geom_text(label = W$sign, 
                y = max(W$value) + 0.1, 
                size = 10, 
                show.legend = F)
  }
  
  p <- p +
    ylim(0, max(W$value) + 2) + 
    geom_text(label = W$impn, 
              y = max(W$value) + 0.2, 
              size = 3, 
              show.legend = F) +
    geom_text(label = W$impper,
              y = max(W$value) + 0.4,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impn_CTRL,
              y = max(W$value) + 0.6,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impper_CTRL,
              y = max(W$value) + 0.8,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impn_MSA,
              y = max(W$value) + 1,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impper_MSA,
              y = max(W$value) + 1.2,
              size = 3,
              show.legend = F) +
  geom_text(label = W$impn_PD,
              y = max(W$value) + 1.4,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impper_PD,
              y = max(W$value) + 1.6,
              size = 3,
              show.legend = F) +
  geom_text(label = W$impn_PSP,
              y = max(W$value) + 1.8,
              size = 3,
              show.legend = F) +
    geom_text(label = W$impper_PSP,
              y = max(W$value) + 2,
              size = 3,
              show.legend = F)
  return(p)
}

Sample-wise contribution to views (areas)

Function

pcs <- function (object, samples = "all", group_by = NULL, return_data = FALSE, ...) 
{
  if (!is(object, "MOFA")) 
    stop("'object' has to be an instance of MOFA")
  scores <- object@cache$contribution_scores %>% 
    reshape2::melt() 
  colnames(scores) <- c("sample", "view", "value")
  
  if (is.null(group_by)) {
    to.plot <- scores
    if (return_data) 
      return(to.plot)
    p <- ggplot(to.plot, aes(x = .data$view, y = .data$value)) + 
      geom_bar(aes(fill = view), stat = "identity", color = "black") + 
      facet_wrap(~sample) + labs(x = "", y = "Contribution score") + 
      theme_classic() + theme(axis.text.x = element_blank(), 
                              axis.ticks.x = element_blank(), legend.position = "top", 
                              legend.title = element_blank())
    return(p)
  }
  else {
    to.plot <- merge(scores, object@samples_metadata[, c("sample", 
                                                         group_by)], by = "sample")
    if (return_data) 
      return(to.plot)
    p <- ggplot(to.plot, aes(x = .data$view, y = .data$value)) + 
      geom_boxplot(aes(fill = view)) + facet_wrap(as.formula(paste("~", 
                                                                   group_by))) + labs(x = "", y = "Contribution score") + 
      theme_classic() + theme(axis.text.x = element_blank(), 
                              axis.ticks.x = element_blank(), legend.position = "top", 
                              legend.title = element_blank())
    return(p)
  }
}

Region-based, with NA, with imputation

https://raw.githack.com/bioFAM/MOFA2_tutorials/master/R_tutorials/downstream_analysis.html

Load model

MOFAobject.trained <- load_model("/work/proteomics/02_data/mofa_nogroups_withna_withimp.hdf5")
## The value -2^31 was detected in the dataset.
## This has been converted to NA within R.
## Warning in .quality_control(object, verbose = verbose): Factor(s) 1, 3 are strongly correlated with the total number of expressed features for at least one of your omics. Such factors appear when there are differences in the total 'levels' between your samples, *sometimes* because of poor normalisation in the preprocessing steps.

Factors correlations

plot_factor_cor(MOFAobject.trained)

Explained variance per area

plot_variance_explained(MOFAobject.trained)

MOFAobject.trained@cache$variance_explained$r2_per_factor
## $group1
##                   CER         FC         MED           SN          STR
## Factor1   2.339857817 2.41563320 2.753508091  2.576148510 24.133199453
## Factor2   0.067317486 5.63617349 7.398545742  7.600158453  0.005507469
## Factor3   0.015407801 0.23275614 0.007551908 13.743638992  0.852268934
## Factor4  11.814665794 0.04336834 0.039535761  0.004839897  2.210634947
## Factor5   4.032844305 2.28035450 0.046372414  2.125096321  5.322384834
## Factor6   0.014066696 2.31698751 8.772975206  1.247757673  0.294625759
## Factor7   1.019167900 0.17645359 1.351934671  5.049353838  1.239514351
## Factor8   2.159816027 0.01425147 2.627253532  1.513928175  2.204763889
## Factor9   0.418668985 1.20673180 0.087374449  0.857770443  4.263985157
## Factor10  0.006401539 0.75365901 0.007772446  0.007688999  4.487466812
plot_variance_explained(MOFAobject.trained,
                        plot_total = T)[[2]]

calculate_variance_explained_per_sample(MOFAobject.trained)$group1 %>%
  as.data.frame() %>% 
  mutate(diagnosis = MOFAobject.trained@samples_metadata$Diagnosis) %>% 
  reshape2::melt(id.vars = "diagnosis") %>% 
  ggplot(aes(variable, value, fill = diagnosis)) +
  geom_boxplot() + 
  geom_point(position = position_jitterdodge(jitter.width = 0.05, dodge.width = 0.7)) +
  theme_bw() + 
  theme(line = element_blank()) +
  labs(x = "", y = "Variance explained", fill = "Diagnosis") +
  scale_fill_manual(values = ggsci::pal_npg()(4))
## Warning: Removed 27 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).

Covariate association

correlate_factors_with_covariates(MOFAobject.trained, 
                                  covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
                                  plot="r",
                                  return_data = F)
## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...

correlate_factors_with_covariates(MOFAobject.trained, 
                                  covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"), 
                                  return_data = F, 
                                  alpha = 1)
## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...

Factor values

plot_factor(MOFAobject.trained, 
            factor = 1:10,
            color_by = "Subtype",
            shape_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factor = 1:10,
            color_by = "Subtype",
            shape_by = "Diagnosis", 
            group_by = "Diagnosis", 
            dodge = T, 
            add_boxplot = T, 
)

Factor 1

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Age",
            shape_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "PMI",
            shape_by = "Diagnosis",
            dot_size = 3
)

Factor 2

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "PMI",
            shape_by = "Diagnosis",
            dot_size = 3,
)

Factor 3

plot_factor(MOFAobject.trained, 
            factors = 3,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

Factor 4

plot_factor(MOFAobject.trained, 
            factors = 4,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

Factor 5

plot_factor(MOFAobject.trained, 
            factors = 5,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 5,
            color_by = "Diagnosis",
            dot_size = 3
)

Factor 7

plot_factor(MOFAobject.trained, 
            factors = 7,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 7,
            color_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 7,
            color_by = "Sex",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 7,
            color_by = "Age",
            dot_size = 3
)

Combined

plot_factors(MOFAobject.trained, 
             factors = 1:3,
             color_by = "Diagnosis",
             shape_by = "Subtype"
)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## Warning: Groups with fewer than two data points have been dropped.
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning: Groups with fewer than two data points have been dropped.
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning: Groups with fewer than two data points have been dropped.
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

plot_factors(MOFAobject.trained, 
             factors = c(1,3),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey")

plot_factors(MOFAobject.trained, 
             factors = c(1,2),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey")

# geom_hline(yintercept = 1, linetype = "dashed", color = "grey")
plot_factors(MOFAobject.trained, 
             factors = c(2,4),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey")

Weights

se.imp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$dat.clean
se.noimp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$dat.norm
meta <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")$expDesign

areas <- names(se.imp@elementMetadata@listData)[grep("imputed", names(se.imp@elementMetadata@listData))] %>% 
  strsplit("_") %>% 
  sget(1)

imp.logi <- se.imp@elementMetadata@listData %>% 
  .[grepl("imputed", names(.))] %>% 
  setNames(areas)

F1, STR

plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
                     view = "STR",
                     factor = 1,
                     nfeatures = 30, 
                     scale = F
)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_data_heatmap(MOFAobject.trained,
                  view = "STR",         # view of interest
                  factor = 1,             # factor of interest
                  features = 15,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = F, 
                  annotation_samples = "Subtype",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "STARD10_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "UTF1_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "ELK4_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "STR",
                  factor = 1, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

F2

Shared proteins

comp <- c("FC", "MED", "SN") %>% 
  lapply(\(aa) {
    tmp <- MOFAobject.trained@expectations$W[[aa]] %>% 
      as.data.frame() %>% 
      dplyr::select(Factor2) %>% 
      arrange(-abs(Factor2))
    
    out <- list()
    
    for (i in c(20, 50, 100)) {
      out[[paste(aa, i, sep = "_")]] <- dplyr::slice(tmp, seq(i)) %>% 
        rownames() %>% 
        strsplit("_") %>% 
        sget(1)
    }
    
    return(out)
  })

seq(3) %>% 
  lapply(\(x) lget(comp, x) %>% 
           setNames(c("FC", "MED", "SN")) %>% 
           fromList() %>% 
           upset(nsets = 3, title = paste0("Top ", c(20, 50, 100)[x])))
## [[1]]

## 
## [[2]]

## 
## [[3]]

lget(comp, 3) %>% 
  setNames(c("FC", "MED", "SN")) %>% 
  unlist() %>% 
  table() %>% 
  data.frame() %>% 
  setNames(c("symbol", "Freq")) %>% 
  filter(Freq >= 2) %>% 
  arrange(-Freq, symbol)
##     symbol Freq
## 1    LIMS1    3
## 2   LRPPRC    3
## 3   SUCLA2    3
## 4   SUCLG1    3
## 5    ACOT8    2
## 6    ALG13    2
## 7     CHN1    2
## 8   COMMD2    2
## 9     CTSG    2
## 10     DCK    2
## 11    FSD1    2
## 12  GDPGP1    2
## 13   GLRX5    2
## 14     GLS    2
## 15    H2AX    2
## 16    H3-7    2
## 17    H4C4    2
## 18   HARS2    2
## 19    HLCS    2
## 20   IARS2    2
## 21  IMPACT    2
## 22   LARS2    2
## 23   LMCD1    2
## 24   MCCC1    2
## 25    NCK1    2
## 26    NCK2    2
## 27 NCKAP1L    2
## 28   NSUN2    2
## 29   NUDT4    2
## 30   NUDT9    2
## 31   PDE5A    2
## 32   PNPT1    2
## 33     POR    2
## 34   PTGR3    2
## 35    RELA    2
## 36  RPRD1B    2
## 37   SMAP2    2
## 38 STARD10    2
## 39  WDSUB1    2

FC

plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
                     view = "FC",
                     factor = 2,
                     nfeatures = 30, 
                     scale = F
)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_data_heatmap(MOFAobject.trained,
                  view = "FC",         # view of interest
                  factor = 2,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained,
            factors = 2, 
            color_by = "LRPPRC_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "SUCLA2_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "PRTN3_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "FC",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

MED

plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
                     view = "MED",
                     factor = 2,
                     nfeatures = 30, 
                     scale = F
)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_data_heatmap(MOFAobject.trained,
                  view = "MED",         # view of interest
                  factor = 2,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "TMEM119_MED", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "SUCLA2_MED", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "MED",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

SN

plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
                     view = "SN",
                     factor = 2,
                     nfeatures = 82, 
                     scale = F
)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_data_heatmap(MOFAobject.trained,
                  view = "SN",         # view of interest
                  factor = 2,             # factor of interest
                  features = 30,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "STARD10_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "RPL34_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "SN",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

F3, SN

plot_top_weights_imp(MOFAobject.trained, imp.list = imp.logi,
                     view = "SN",
                     factor = 3,
                     nfeatures = 30, 
                     scale = F
)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

plot_data_heatmap(MOFAobject.trained,
                  view = "SN",         # view of interest
                  factor = 3,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 3, 
            color_by = "RPL32_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 3, 
            color_by = "PDE1B_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "SN",
                  factor = 3, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

tSNE

model <- run_tsne(MOFAobject.trained, perplexity = 10)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Subtype", 
            shape_by = "Diagnosis", 
            dot_size = 3
)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the MOFA2 package.
##   Please report the issue at <https://github.com/bioFAM/MOFA2>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Diagnosis", 
            # shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor1", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor2", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor3", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

Sample-wise contribution scores

MOFAobject.trained <- calculate_contribution_scores(object = MOFAobject.trained)
samples_metadata(MOFAobject.trained)
##    Age Diagnosis Duration Hemisphere Origin   PMI    Sex  Subtype  group sample
## 1   65       PSP        7       Left    BBH  28.0   Male       RS group1   1370
## 2   73       MSA        7       Left    BBH  24.0   Male        P group1   1379
## 3   74       MSA        4       Left    BBH  56.0 Female        C group1   1388
## 4   59       MSA        8      Right    BBH  24.0 Female        C group1   1391
## 5   56       MSA        8      Right    BBH  72.0   Male      C+P group1   1392
## 6   80       PSP        9      Right    BBH  36.0   Male       RS group1   1394
## 7   55       MSA        3       Left    BBH  31.0 Female        C group1   1398
## 8   66        PD       13      Right    BBH  53.0   Male     None group1   1400
## 9   69       MSA        8       Left    BBH  24.0 Female        P group1   1406
## 10  56       MSA        7      Right    BBH  58.0   Male        P group1   1407
## 11  70       MSA        9      Right    BBH  26.0   Male        P group1   1436
## 12  79        PD       10       Left    BBH 118.0   Male     None group1   1446
## 13  68       MSA        7      Right    BBH  70.0   Male        P group1   1462
## 14  91      CTRL       NA       Left    BBH  24.0 Female     None group1   1464
## 15  68      CTRL       NA      Right    BBH  24.0   Male     None group1   1467
## 16  91      CTRL       NA       Left    BBH  24.0 Female     None group1   1490
## 17  82      CTRL       NA      Right    BBH  48.0   Male     None group1   1492
## 18  86      CTRL       NA       Left    BBH  45.0 Female     None group1   1494
## 19  97      CTRL       NA      Right    BBH  90.0 Female     None group1   1495
## 20  82       PSP        5       Left    BBH  24.0   Male     None group1   1385
## 21  87       PSP        8      Right    BBH  30.0 Female     None group1   1387
## 22  67       PSP        9       Left    BBH  30.0   Male     None group1   1397
## 23  71       PSP        4      Right    BBH  30.0   Male       RS group1   1401
## 24  69       PSP        9      Right    BBH  29.0   Male Cortical group1   1410
## 25  81        PD        9       Left    BBH  34.0   Male     None group1   1411
## 26  70       PSP        6       Left    BBH  19.0   Male       RS group1   1414
## 27  71       PSP        8      Right    BBH 101.0   Male Cortical group1   1420
## 28  64        PD       13       Left    BBH 117.5   Male     None group1   1440
## 29  88      CTRL       NA      Right    BBH  24.0 Female     None group1   1453
## 30  72        PD        6      Right    BBH  89.0   Male     None group1   1457
## 31  84      CTRL       NA      Right    BBH  96.0 Female     None group1   1466
## 32  75       PSP        7      Right    BBH  15.0 Female       RS group1   1381
## 33  95        PD        7       Left    BBH  86.0 Female     None group1   1412
## 34  91        PD       10       Left    BBH  39.0 Female     None group1   1458
## 35  86        PD       12       Left    BBH  39.0   Male     None group1   1424
## 36  75        PD       NA       Left    BBH  48.0   Male     None group1   1449
## 37  48      CTRL       NA      Whole    BBH  48.0 Female     None group1   1461
## 38  91      CTRL       NA       Left    BBH  96.0 Female     None group1   1465
## 39  75        PD       NA       Left    BBH  32.0 Female     None group1   1429
## 40  62        PD       NA      Right    BBH  34.0   Male     None group1   1469
##    CER_contribution FC_contribution MED_contribution SN_contribution
## 1        0.15082781      0.06216384       0.12712092      0.46087914
## 2        0.14714621      0.27634369       0.17855782      0.17427280
## 3        0.12335640      0.22708067       0.21952391      0.29500351
## 4               NaN             NaN              NaN             NaN
## 5               NaN             NaN              NaN             NaN
## 6        0.22138774      0.13635375       0.23508441      0.24971132
## 7               NaN             NaN              NaN             NaN
## 8        0.23427308      0.24046448       0.20627061      0.16327756
## 9        0.23135490      0.21387768       0.14590016      0.16719717
## 10              NaN             NaN              NaN             NaN
## 11              NaN             NaN              NaN             NaN
## 12       0.21281771      0.23732481       0.14386542      0.22032829
## 13              NaN             NaN              NaN             NaN
## 14       0.26906764      0.25865746       0.16421821      0.13261004
## 15              NaN             NaN              NaN             NaN
## 16              NaN             NaN              NaN             NaN
## 17       0.13177486      0.16031751       0.03116276      0.49728780
## 18       0.22941892      0.22196431       0.05917558      0.31262607
## 19       0.18399271      0.20715479       0.27007271      0.19715939
## 20              NaN             NaN              NaN             NaN
## 21       0.25899102      0.00000000       0.14081864      0.16015999
## 22       0.32527763      0.11095035       0.08469555      0.27869658
## 23       0.26404507      0.10151794       0.10497410      0.33232702
## 24              NaN             NaN              NaN             NaN
## 25       0.09431939      0.01175833       0.10577559      0.63934940
## 26       0.27749772      0.07444969       0.26094747      0.13652442
## 27       0.37322884      0.00000000       0.23567577      0.14646064
## 28       0.10797343      0.32450990       0.27192454      0.23417692
## 29       0.16059805      0.08809670       0.28351085      0.24027018
## 30       0.21043819      0.32227562       0.15072773      0.15046156
## 31       0.35634937      0.15094665       0.22551296      0.09514165
## 32              NaN             NaN              NaN             NaN
## 33              NaN             NaN              NaN             NaN
## 34              NaN             NaN              NaN             NaN
## 35              NaN             NaN              NaN             NaN
## 36              NaN             NaN              NaN             NaN
## 37              NaN             NaN              NaN             NaN
## 38              NaN             NaN              NaN             NaN
## 39              NaN             NaN              NaN             NaN
## 40              NaN             NaN              NaN             NaN
##    STR_contribution
## 1        0.19900829
## 2        0.22367948
## 3        0.13503551
## 4               NaN
## 5               NaN
## 6        0.15746278
## 7               NaN
## 8        0.15571427
## 9        0.24167008
## 10              NaN
## 11              NaN
## 12       0.18566378
## 13              NaN
## 14       0.17544665
## 15              NaN
## 16              NaN
## 17       0.17945707
## 18       0.17681512
## 19       0.14162040
## 20              NaN
## 21       0.44003035
## 22       0.20037988
## 23       0.19713586
## 24              NaN
## 25       0.14879729
## 26       0.25058070
## 27       0.24463476
## 28       0.06141521
## 29       0.22752422
## 30       0.16609690
## 31       0.17204937
## 32              NaN
## 33              NaN
## 34              NaN
## 35              NaN
## 36              NaN
## 37              NaN
## 38              NaN
## 39              NaN
## 40              NaN
pcs(MOFAobject.trained, group_by = "group")
## Warning: Removed 95 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

pcs(MOFAobject.trained)
## Warning: Removed 95 rows containing missing values or values outside the scale range
## (`geom_bar()`).

Region-based, with NA, no imputation

https://raw.githack.com/bioFAM/MOFA2_tutorials/master/R_tutorials/downstream_analysis.html

Load model

MOFAobject.trained <- load_model("/work/proteomics/02_data/mofa_nogroups_withna_noimp.hdf5")
## The value -2^31 was detected in the dataset.
## This has been converted to NA within R.
## Warning in .quality_control(object, verbose = verbose): Factor(s) 1, 3 are strongly correlated with the total number of expressed features for at least one of your omics. Such factors appear when there are differences in the total 'levels' between your samples, *sometimes* because of poor normalisation in the preprocessing steps.

Factors correlations

plot_factor_cor(MOFAobject.trained)

# =========================
# Setup
# =========================
W_list <- MOFAobject.trained@expectations$W  # list of view matrices (features x factors)

# Fixed view order (change if you prefer a different order)
view_order <- c("CER", "FC", "MED", "SN", "STR")
W_list <- W_list[view_order]

# Ensure factor names exist and are consistent (Factor1..FactorN)
normalize_factor_names <- function(m) {
  if (is.null(colnames(m)) || any(nchar(colnames(m)) == 0)) {
    colnames(m) <- paste0("Factor", seq_len(ncol(m)))
  }
  m
}
W_list <- lapply(W_list, normalize_factor_names)

# Helper: strip trailing "_VIEW" suffix (e.g., IGKV3-7_CER -> IGKV3-7)
strip_suffix <- function(ids) sub("_[^_]+$", "", ids)

# Pairwise correlations of W (features x factors) between two views on common features
pair_cor_views <- function(W_i, W_j, method = "pearson") {
  rid_i <- strip_suffix(rownames(W_i))
  rid_j <- strip_suffix(rownames(W_j))
  common <- intersect(rid_i, rid_j)

  # If overlap is tiny, return NA matrix with proper dimnames
  if (length(common) < 3) {
    out <- matrix(NA_real_, nrow = ncol(W_i), ncol = ncol(W_j))
    rownames(out) <- colnames(W_i)
    colnames(out) <- colnames(W_j)
    return(out)
  }

  A <- W_i[match(common, rid_i), , drop = FALSE]
  B <- W_j[match(common, rid_j), , drop = FALSE]

  out <- cor(A, B, use = "pairwise.complete.obs", method = method)
  # Preserve factor names
  rownames(out) <- colnames(W_i)
  colnames(out) <- colnames(W_j)
  out
}

# =========================
# Compute only upper triangle (including diagonal) panels
# =========================
method  <- "pearson"  # or "spearman"
top_k   <- 3          # number of top correlations to annotate per panel

pairwise_mats <- list()
for (i in seq_along(view_order)) {
  for (j in seq_along(view_order)) {
    if (i <= j) {  # upper triangle + diagonal
      v1 <- view_order[i]; v2 <- view_order[j]
      lbl <- paste(v1, v2, sep = " vs ")
      pairwise_mats[[lbl]] <- pair_cor_views(W_list[[v1]], W_list[[v2]], method = method)
    }
  }
}

# =========================
# Melt for ggplot facets
# =========================
if (!requireNamespace("reshape2", quietly = TRUE)) install.packages("reshape2")
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")

melt_pair <- function(mat, v1, v2) {
  df <- reshape2::melt(mat)
  colnames(df) <- c("Factor1", "Factor2", "cor")
  df$view1 <- v1
  df$view2 <- v2
  df
}

df_all <- do.call(
  rbind,
  lapply(names(pairwise_mats), function(lbl) {
    parts <- strsplit(lbl, " vs ", fixed = TRUE)[[1]]
    melt_pair(pairwise_mats[[lbl]], parts[1], parts[2])
  })
)

# Factor ordering (ensures axes read Factor1..FactorN everywhere)
all_factors <- sort(unique(unlist(lapply(W_list, colnames))))
df_all$Factor1 <- factor(df_all$Factor1, levels = all_factors)
df_all$Factor2 <- factor(df_all$Factor2, levels = all_factors)

# =========================
# Build top-k overlay (per panel) EXCLUDING r == 1 and diagonal self-pairs
# =========================
split_panels <- split(df_all, interaction(df_all$view1, df_all$view2, drop = TRUE))

top_list <- lapply(split_panels, function(d) {
  # Drop NA and exact 1.0 correlations
  d <- subset(d, !is.na(cor) & cor != 1)

  # In within-view panels, exclude diagonal cells (Factor1 == Factor2)
  # which would be 1 by definition; already excluded by cor != 1,
  # but keep the explicit filter in case of numerical jitter.
  if (nrow(d) > 0 && unique(d$view1) == unique(d$view2)) {
    d <- subset(d, as.character(Factor1) != as.character(Factor2))
  }

  # Order by absolute correlation descending
  d <- d[order(-abs(d$cor)), , drop = FALSE]

  # Keep top-k
  if (nrow(d) == 0) return(d)
  d[seq_len(min(nrow(d), top_k)), , drop = FALSE]
})

top_df <- do.call(rbind, top_list)

# =========================
# Plot
# =========================
labeller_fn <- ggplot2::labeller(
  view1 = function(x) paste0("View 1: ", x),
  view2 = function(x) paste0("View 2: ", x)
)

p <- ggplot2::ggplot(df_all, ggplot2::aes(x = Factor2, y = Factor1, fill = cor)) +
  ggplot2::geom_tile() +
  # Overlay: mark top-k tiles with a ring + correlation value
  ggplot2::geom_point(
    data = top_df,
    shape = 21, stroke = 0.9, size = 3.8, colour = "black", fill = NA
  ) +
  ggplot2::geom_text(
    data = top_df,
    ggplot2::aes(label = sprintf("%.2f", cor)),
    size = 3, fontface = "bold", colour = "black"
  ) +
  ggplot2::scale_fill_gradient2(
    low = "navy", mid = "white", high = "firebrick3",
    midpoint = 0, limits = c(-1, 1), oob = scales::squish, name = "r"
  ) +
  ggplot2::facet_grid(view1 ~ view2, labeller = labeller_fn, drop = TRUE) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(face = "bold"),
    axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
    panel.spacing = grid::unit(0.6, "lines"),
    legend.position = "right"
  ) +
  ggplot2::labs(
    title    = sprintf("MOFA W: Factor × Factor correlations (upper triangle), method = %s", method),
    subtitle = sprintf("Top-%d factor pairs per panel highlighted (excluding r = 1); panels use pairwise intersections of features", top_k),
    x = "Factors in View 2",
    y = "Factors in View 1"
  )

print(p)

Explained variance per area

plot_variance_explained(MOFAobject.trained)

MOFAobject.trained@cache$variance_explained$r2_per_factor
## $group1
##                  CER          FC          MED          SN          STR
## Factor1  0.072622299 3.231221437  2.674627304  0.01284480 36.455994844
## Factor2 14.094030857 3.056883812  9.029161930  1.87025070  5.805820227
## Factor3  6.209421158 3.535330296  0.028491020 19.61889863  1.466298103
## Factor4  6.710535288 5.907905102  0.435268879  2.19438076 10.662490129
## Factor5  0.005841255 1.038879156 17.748779058  0.00308156  4.222744703
## Factor6  2.574747801 1.141369343  6.447696686 11.38798594  0.010675192
## Factor7  7.396245003 6.605488062  0.034427643  7.38848448  0.007265806
## Factor8  0.006121397 0.005877018  0.008529425 10.15772820  0.357550383
plot_variance_explained(MOFAobject.trained,
                        plot_total = T)[[2]]

Covariate association

correlate_factors_with_covariates(MOFAobject.trained, 
                                  covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"),
                                  plot="r",
                                  return_data = F)
## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...

correlate_factors_with_covariates(MOFAobject.trained, 
                                  covariates = c("Diagnosis", "Duration", "Hemisphere", "PMI", "Sex", "Subtype", "Age"), 
                                  return_data = F, 
                                  alpha = 1)
## Warning in correlate_factors_with_covariates(MOFAobject.trained, covariates =
## c("Diagnosis", : There are non-numeric values in the covariates data.frame,
## converting to numeric...

Factor values

plot_factor(MOFAobject.trained, 
            factor = 1:8,
            color_by = "Subtype",
            shape_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factor = 1:8,
            color_by = "Subtype",
            shape_by = "Diagnosis", 
            group_by = "Diagnosis", 
            dodge = T, 
            add_boxplot = T, 
)

Factor 1

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "Age",
            shape_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 1,
            color_by = "PMI",
            shape_by = "Diagnosis",
            dot_size = 3
)

Factor 2

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "Diagnosis",
            dot_size = 3
)

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

plot_factor(MOFAobject.trained, 
            factors = 2,
            color_by = "Sex",
            shape_by = "Diagnosis",
            dot_size = 3,
)

Factor 3

plot_factor(MOFAobject.trained, 
            factors = 3,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

Factor 4

plot_factor(MOFAobject.trained, 
            factors = 4,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

Factor 5

plot_factor(MOFAobject.trained, 
            factors = 5,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis"
)

Factor 6

plot_factor(MOFAobject.trained, 
            factors = 6,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis", 
)

plot_factor(MOFAobject.trained, 
            factors = 6,
            color_by = "Diagnosis",
            dot_size = 3
)

Factor 8

plot_factor(MOFAobject.trained, 
            factors = 8,
            color_by = "Subtype",
            shape_by = "Diagnosis",
            dot_size = 3,
            group_by = "Diagnosis", 
)

plot_factor(MOFAobject.trained, 
            factors = 8,
            color_by = "Diagnosis",
            dot_size = 3
)

Combined

plot_factors(MOFAobject.trained, 
             factors = c(1,3),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 0.9, linetype = "dashed", color = "grey") +
  geom_hline(yintercept = 0.9, linetype = "dashed", color = "grey")

plot_factors(MOFAobject.trained, 
             factors = c(1,2),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey") +
  geom_hline(yintercept = -1.1, linetype = "dashed", color = "grey")

plot_factors(MOFAobject.trained, 
             factors = c(6,8),
             color_by = "Subtype",
             shape_by = "Diagnosis",
             dot_size = 3) +
  geom_vline(xintercept = 1.15, linetype = "dashed", color = "grey")

# geom_hline(yintercept = 1, linetype = "dashed", color = "grey")

Weights

F1, STR

plot_top_weights(MOFAobject.trained,
                 view = "STR",
                 factor = 1,
                 nfeatures = 30, 
                 scale = F
)

plot_data_heatmap(MOFAobject.trained,
                  view = "STR",         # view of interest
                  factor = 1,             # factor of interest
                  features = 15,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = F, 
                  annotation_samples = "Subtype",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "PDE1B_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "UTF1_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 1, 
            color_by = "ATCAY_STR", 
            group_by = "Subtype",
            shape_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "STR",
                  factor = 1, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

F2

Shared proteins

comp <- c("FC", "MED", "SN") %>% 
  lapply(\(aa) {
    tmp <- MOFAobject.trained@expectations$W[[aa]] %>% 
      as.data.frame() %>% 
      dplyr::select(Factor2) %>% 
      arrange(-abs(Factor2))
    
    out <- list()
    
    for (i in c(20, 50, 100)) {
      out[[paste(aa, i, sep = "_")]] <- dplyr::slice(tmp, seq(i)) %>% 
        rownames() %>% 
        strsplit("_") %>% 
        sget(1)
    }
    
    return(out)
  })

seq(3) %>% 
  lapply(\(x) lget(comp, x) %>% 
           setNames(c("FC", "MED", "SN")) %>% 
           fromList() %>% 
           upset(nsets = 3, title = paste0("Top ", c(20, 50, 100)[x])))
## [[1]]

## 
## [[2]]

## 
## [[3]]

FC

plot_top_weights(MOFAobject.trained,
                 view = "FC",
                 factor = 2,
                 nfeatures = 30, 
                 scale = F
)

plot_data_heatmap(MOFAobject.trained,
                  view = "FC",         # view of interest
                  factor = 2,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "PTGR3_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "OGDHL_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "GRPEL1_FC", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "FC",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

MED

plot_top_weights(MOFAobject.trained,
                 view = "MED",
                 factor = 2,
                 nfeatures = 30, 
                 scale = F
)

plot_data_heatmap(MOFAobject.trained,
                  view = "MED",         # view of interest
                  factor = 2,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "FSD1_MED", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "SUCLA2_MED", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "MED",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

SN

plot_top_weights(MOFAobject.trained,
                 view = "SN",
                 factor = 2,
                 nfeatures = 30, 
                 scale = F
)

plot_data_heatmap(MOFAobject.trained,
                  view = "SN",         # view of interest
                  factor = 2,             # factor of interest
                  features = 30,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "CBR1_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 2, 
            color_by = "TRAP1_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "SN",
                  factor = 2, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

F3, SN

plot_top_weights(MOFAobject.trained,
                 view = "SN",
                 factor = 3,
                 nfeatures = 30, 
                 scale = F
)

plot_data_heatmap(MOFAobject.trained,
                  view = "SN",         # view of interest
                  factor = 3,             # factor of interest
                  features = 20,          # number of features to plot (they are selected by weight)
                  cluster_rows = T, 
                  cluster_cols = T,
                  show_rownames = TRUE, 
                  show_colnames = T, 
                  annotation_samples = "Diagnosis",
                  denoise = T
)

plot_factor(MOFAobject.trained, 
            factors = 3, 
            color_by = "XPO1_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_factor(MOFAobject.trained, 
            factors = 3, 
            color_by = "CSE1L_SN", 
            group_by = "Diagnosis",
            add_violin = TRUE,
)

plot_data_scatter(MOFAobject.trained, 
                  view = "SN",
                  factor = 3, 
                  features = 4,
                  color_by = "Subtype",
                  shape_by = "Diagnosis"
)

tSNE

model <- run_tsne(MOFAobject.trained, perplexity = 10)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Subtype", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Diagnosis", 
            # shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor1", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor2", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

plot_dimred(model,
            method = "TSNE",  # method can be either "TSNE" or "UMAP"
            color_by = "Factor3", 
            shape_by = "Diagnosis", 
            dot_size = 3
)

Sample-wise contribution scores

MOFAobject.trained <- calculate_contribution_scores(object = MOFAobject.trained)
samples_metadata(MOFAobject.trained)
##    Age Diagnosis Duration Hemisphere Origin   PMI    Sex  Subtype  group sample
## 1   65       PSP        7       Left    BBH  28.0   Male       RS group1   1370
## 2   73       MSA        7       Left    BBH  24.0   Male        P group1   1379
## 3   74       MSA        4       Left    BBH  56.0 Female        C group1   1388
## 4   59       MSA        8      Right    BBH  24.0 Female        C group1   1391
## 5   56       MSA        8      Right    BBH  72.0   Male      C+P group1   1392
## 6   80       PSP        9      Right    BBH  36.0   Male       RS group1   1394
## 7   55       MSA        3       Left    BBH  31.0 Female        C group1   1398
## 8   66        PD       13      Right    BBH  53.0   Male     None group1   1400
## 9   69       MSA        8       Left    BBH  24.0 Female        P group1   1406
## 10  56       MSA        7      Right    BBH  58.0   Male        P group1   1407
## 11  70       MSA        9      Right    BBH  26.0   Male        P group1   1436
## 12  79        PD       10       Left    BBH 118.0   Male     None group1   1446
## 13  68       MSA        7      Right    BBH  70.0   Male        P group1   1462
## 14  91      CTRL       NA       Left    BBH  24.0 Female     None group1   1464
## 15  68      CTRL       NA      Right    BBH  24.0   Male     None group1   1467
## 16  91      CTRL       NA       Left    BBH  24.0 Female     None group1   1490
## 17  82      CTRL       NA      Right    BBH  48.0   Male     None group1   1492
## 18  86      CTRL       NA       Left    BBH  45.0 Female     None group1   1494
## 19  97      CTRL       NA      Right    BBH  90.0 Female     None group1   1495
## 20  82       PSP        5       Left    BBH  24.0   Male     None group1   1385
## 21  87       PSP        8      Right    BBH  30.0 Female     None group1   1387
## 22  67       PSP        9       Left    BBH  30.0   Male     None group1   1397
## 23  71       PSP        4      Right    BBH  30.0   Male       RS group1   1401
## 24  69       PSP        9      Right    BBH  29.0   Male Cortical group1   1410
## 25  81        PD        9       Left    BBH  34.0   Male     None group1   1411
## 26  70       PSP        6       Left    BBH  19.0   Male       RS group1   1414
## 27  71       PSP        8      Right    BBH 101.0   Male Cortical group1   1420
## 28  64        PD       13       Left    BBH 117.5   Male     None group1   1440
## 29  88      CTRL       NA      Right    BBH  24.0 Female     None group1   1453
## 30  72        PD        6      Right    BBH  89.0   Male     None group1   1457
## 31  84      CTRL       NA      Right    BBH  96.0 Female     None group1   1466
## 32  75       PSP        7      Right    BBH  15.0 Female       RS group1   1381
## 33  95        PD        7       Left    BBH  86.0 Female     None group1   1412
## 34  91        PD       10       Left    BBH  39.0 Female     None group1   1458
## 35  86        PD       12       Left    BBH  39.0   Male     None group1   1424
## 36  75        PD       NA       Left    BBH  48.0   Male     None group1   1449
## 37  48      CTRL       NA      Whole    BBH  48.0 Female     None group1   1461
## 38  91      CTRL       NA       Left    BBH  96.0 Female     None group1   1465
## 39  75        PD       NA       Left    BBH  32.0 Female     None group1   1429
## 40  62        PD       NA      Right    BBH  34.0   Male     None group1   1469
##    CER_contribution FC_contribution MED_contribution SN_contribution
## 1        0.14214162      0.07123318        0.2335368       0.3917811
## 2        0.06725632      0.31186582        0.2231970       0.1263306
## 3        0.18876339      0.21377168        0.2352642       0.2225193
## 4               NaN             NaN              NaN             NaN
## 5               NaN             NaN              NaN             NaN
## 6        0.16822988      0.23762093        0.1425290       0.2873132
## 7               NaN             NaN              NaN             NaN
## 8        0.22864984      0.23697852        0.2687660       0.1936276
## 9        0.14957799      0.11660456        0.2413496       0.1457976
## 10              NaN             NaN              NaN             NaN
## 11              NaN             NaN              NaN             NaN
## 12       0.21213553      0.26038667        0.1666308       0.1829730
## 13              NaN             NaN              NaN             NaN
## 14       0.26612438      0.20493846        0.2456342       0.1432734
## 15              NaN             NaN              NaN             NaN
## 16              NaN             NaN              NaN             NaN
## 17       0.17166178      0.14881793        0.2031213       0.3280382
## 18       0.16583917      0.16992017        0.0645116       0.4507353
## 19       0.18931714      0.27299803        0.1633389       0.1954027
## 20              NaN             NaN              NaN             NaN
## 21       0.25027903      0.07958009        0.1626316       0.3382785
## 22       0.29430106      0.04239646        0.1040508       0.3707292
## 23       0.13338611      0.16168388        0.1905265       0.2708038
## 24              NaN             NaN              NaN             NaN
## 25       0.18935370      0.00000000        0.1647600       0.5344074
## 26       0.27797464      0.00000000        0.2874644       0.2305561
## 27       0.35424787      0.02548541        0.1898440       0.2271896
## 28       0.19559509      0.24869516        0.2153742       0.2067256
## 29       0.15693131      0.05634121        0.3298812       0.2506410
## 30       0.21998408      0.28382531        0.1515171       0.1386265
## 31       0.14490946      0.20805933        0.3463598       0.1873873
## 32              NaN             NaN              NaN             NaN
## 33              NaN             NaN              NaN             NaN
## 34              NaN             NaN              NaN             NaN
## 35              NaN             NaN              NaN             NaN
## 36              NaN             NaN              NaN             NaN
## 37              NaN             NaN              NaN             NaN
## 38              NaN             NaN              NaN             NaN
## 39              NaN             NaN              NaN             NaN
## 40              NaN             NaN              NaN             NaN
##    STR_contribution
## 1        0.16130730
## 2        0.27135021
## 3        0.13968144
## 4               NaN
## 5               NaN
## 6        0.16430709
## 7               NaN
## 8        0.07197804
## 9        0.34667028
## 10              NaN
## 11              NaN
## 12       0.17787409
## 13              NaN
## 14       0.14002964
## 15              NaN
## 16              NaN
## 17       0.14836080
## 18       0.14899375
## 19       0.17894316
## 20              NaN
## 21       0.16923073
## 22       0.18852240
## 23       0.24359972
## 24              NaN
## 25       0.11147885
## 26       0.20400486
## 27       0.20323313
## 28       0.13360999
## 29       0.20620532
## 30       0.20604697
## 31       0.11328414
## 32              NaN
## 33              NaN
## 34              NaN
## 35              NaN
## 36              NaN
## 37              NaN
## 38              NaN
## 39              NaN
## 40              NaN
pcs(MOFAobject.trained, group_by = "group")
## Warning: Removed 95 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

pcs(MOFAobject.trained)
## Warning: Removed 95 rows containing missing values or values outside the scale range
## (`geom_bar()`).

Session info

tt - Sys.time()
## Time difference of -1.257868 mins
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2;  LAPACK version 3.12.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MOFAdata_1.24.0 tibble_3.2.1    purrr_1.0.4     tidyr_1.3.1    
##  [5] qs_0.27.3       scHelper_0.0.5  UpSetR_1.4.0    ggplot2_3.5.2  
##  [9] MOFA2_1.18.0    dplyr_1.1.4     magrittr_2.0.3 
## 
## loaded via a namespace (and not attached):
##   [1] mnormt_2.1.1                gridExtra_2.3              
##   [3] rlang_1.1.6                 clue_0.3-66                
##   [5] GetoptLong_1.0.5            matrixStats_1.5.0          
##   [7] compiler_4.5.2              mgcv_1.9-3                 
##   [9] sccore_1.0.5                dir.expiry_1.16.0          
##  [11] png_0.1-8                   vctrs_0.6.5                
##  [13] reshape2_1.4.4              stringr_1.5.1              
##  [15] pkgconfig_2.0.3             shape_1.4.6.1              
##  [17] crayon_1.5.3                fastmap_1.2.0              
##  [19] backports_1.5.0             XVector_0.48.0             
##  [21] labeling_0.4.3              rmarkdown_2.29             
##  [23] UCSC.utils_1.4.0            xfun_0.52                  
##  [25] cachem_1.1.0                GenomeInfoDb_1.44.0        
##  [27] jsonlite_2.0.0              rhdf5filters_1.20.0        
##  [29] DelayedArray_0.34.1         Rhdf5lib_1.30.0            
##  [31] psych_2.5.3                 broom_1.0.8                
##  [33] parallel_4.5.2              cluster_2.1.8.1            
##  [35] R6_2.6.1                    bslib_0.9.0                
##  [37] stringi_1.8.7               RColorBrewer_1.1-3         
##  [39] reticulate_1.42.0           GGally_2.2.1               
##  [41] car_3.1-3                   leidenAlg_1.1.5            
##  [43] GenomicRanges_1.60.0        jquerylib_0.1.4            
##  [45] Rcpp_1.0.14                 SummarizedExperiment_1.38.1
##  [47] iterators_1.0.14            knitr_1.50                 
##  [49] IRanges_2.42.0              splines_4.5.2              
##  [51] Matrix_1.7-3                igraph_2.1.4               
##  [53] tidyselect_1.2.1            rstudioapi_0.17.1          
##  [55] dichromat_2.0-0.1           abind_1.4-8                
##  [57] yaml_2.3.10                 stringfish_0.16.0          
##  [59] doParallel_1.0.17           codetools_0.2-20           
##  [61] lattice_0.22-7              plyr_1.8.9                 
##  [63] Biobase_2.68.0              basilisk.utils_1.20.0      
##  [65] withr_3.0.2                 evaluate_1.0.3             
##  [67] Rtsne_0.17                  ggstats_0.9.0              
##  [69] RcppParallel_5.1.10         circlize_0.4.16            
##  [71] ggpubr_0.6.0                pillar_1.10.2              
##  [73] filelock_1.0.3              carData_3.0-5              
##  [75] MatrixGenerics_1.20.0       corrplot_0.95              
##  [77] foreach_1.5.2               stats4_4.5.2               
##  [79] generics_0.1.4              S4Vectors_0.46.0           
##  [81] scales_1.4.0                RApiSerialize_0.1.4        
##  [83] glue_1.8.0                  pheatmap_1.0.12            
##  [85] tools_4.5.2                 ggsignif_0.6.4             
##  [87] forcats_1.0.1               cowplot_1.1.3              
##  [89] rhdf5_2.52.0                grid_4.5.2                 
##  [91] conos_1.5.2                 colorspace_2.1-1           
##  [93] GenomeInfoDbData_1.2.14     nlme_3.1-168               
##  [95] basilisk_1.20.0             HDF5Array_1.36.0           
##  [97] Formula_1.2-5               cli_3.6.5                  
##  [99] S4Arrays_1.8.0              ComplexHeatmap_2.24.0      
## [101] uwot_0.2.3                  gtable_0.3.6               
## [103] rstatix_0.7.2               ggsci_3.2.0                
## [105] sass_0.4.10                 digest_0.6.37              
## [107] BiocGenerics_0.54.0         SparseArray_1.8.0          
## [109] ggrepel_0.9.6               rjson_0.2.23               
## [111] farver_2.1.2                htmltools_0.5.8.1          
## [113] lifecycle_1.0.4             httr_1.4.7                 
## [115] h5mread_1.0.0               GlobalOptions_0.1.2